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Abstract 

Motivated by applications in thermal QCD and cosmology, we elaborate on a general method 
for computing next-to-leading order spectral functions for composite operators at vanishing 
spatial momentum, accounting for real, virtual as well as thermal corrections. As an exam- 
ple, we compute these functions (together with the corresponding imaginary-time correlators 
which can be compared with lattice simulations) for scalar and pseudoscalar densities in pure 
Yang-Mills theory. Our results may turn out to be helpful in non-perturbative estimates of 
the corresponding transport coefficients, which are the bulk viscosity in the scalar channel 
and the rate of anomalous chirality violation in the pseudoscalar channel. We also mention 
links to cosmology, although the most useful results in that context may come from a future 
generalization of our methods to other correlators. 



August 2011 



1. Introduction 

A "spectral function", p, which has a Minkowskian four-momentum, V, as its argument, 
can formally be defined as a Fourier transform of the thermal expectation value ((. . .) T = 
-g>TV [e~P H (. . .)]) of the commutator of a given gauge-invariant local operator: 

p(V) = |^d^e^(l[6(t,x),d(0,0)]) T , (1.1) 

where X = (t,x), the operator is defined in the Heisenberg picture, and V ■ X = ut — p • x. 
All other orderings (retarded, time-ordered, Wightman, etc.) can be expressed in terms of 
the spectral function [1, 2]. Therefore various spectral functions play an important role in 
theoretical analyses of thermal systems. In particular, in the "linear response" regime, the 
production rate of a weakly-coupled particle species from the medium is proportional to the 
spectral function of the operator that the particle in question couples to [1, 2]. 

Because of the fundamental role of spectral functions, they make an appearance in many 
different areas of many-body physics; in the following, selected examples having a connection 
to our present study are outlined in some more detail. 

QCD at high temperatures 

Various "transport coefficients" , which characterize the relaxation of local excesses in con- 
served currents (momentum; number densities; electromagnetic current) to their equilibrium 
values, can be viewed as the "low-energy constants" of an interacting hot plasma (cf. e.g. 
ref. [3]). They can be defined through the limit lim tJ _ 5> o+ Pi^i ^)/ LjJ °f an appropriate spectral 
function. In recent years considerable efforts have been devoted to determining these quanti- 
ties at temperatures around a few hundred MeV, relevant for the current experimental heavy 
ion collision program; an extensive review, with a particular perspective on non-perturbative 
lattice measurements, can be found in ref. [4]. 

Unfortunately, as becomes clear from ref. [4], a reliable non-perturbative determination of 
transport coefficients represents a formidable challenge, even in the limit of high temperatures 
where it can be argued that pure Yang-Mills theory should already yield a qualitatively 
correct description. The problem is that direct lattice measurements concern Euclidean 
correlators, whereas transport coefficients are related to the longest Minkowskian time scales 
characterizing the system. Even though the two situations can in principle be related to 
each other through an analytic continuation [5], a necessary requirement appears to be the 
subtraction of short-distance divergences [6] , and even then the problem remains numerically 
ill-posed (see e.g. ref. [4] and references therein). It therefore seems pertinent to obtain as 
much information as possible by analytic means, in order to carry out a subtraction leaving 
over a simple function containing (mostly) only infrared physics. Thanks to asymptotic 
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freedom, it is possible to address analytically not only the temperature-independent short- 
distance divergences but also, at high enough temperatures, thermal modifications originating 
from the large-frequency range uj ~ ttT, and this is one of the goals of the present work. 

QCD at low temperatures 

At low temperatures, where chiral symmetry is spontaneously broken, numerical simulations 
are particularly demanding; on the other hand, chiral perturbation theory [7], if pursued to 
a sufficient order, may provide for a reasonably accurate description of this regime. Since 
mesons are bosons, the analytic computations encountered are technically not unlike those 
at very high temperatures in pure Yang-Mills theory, and many of the same methods may 
turn out to be useful. This then serves as an additional motivation for developing analytic 
methods. Recent computations of various transport coefficients within chiral perturbation 
theory (or extensions thereof) can be found in refs. [8]-[12]; we envisage that our techniques 
and generalizations thereof may allow to address some of the open issues outlined in ref. [8]. 

Cosmology 

In cosmology, particle production rates are relevant for certain Dark Matter scenarios, as 
well as for leptogenesis computations in which a non-thermal distribution of unstable particles 
serves as an intermediate stage. Typically, one is then interested in the whole non-equilibrium 
spectrum of the particles produced, proportional to a spectral function evaluated at V = 
(E p ,p), where p is the spatial momentum and E p = \/ p 2 + M 2 (cf. e.g. ref. [13]). However, 
if the particles produced are non-relativistic, with a mass M»T, then it can be argued that 
the spectral function at V = (M, 0) already yields some of the information. An example of a 
recent analysis in this spirit, in which a subset of the 2-loop topologies studied in the present 
paper were computed to leading non-trivial order in an expansion in T 2 /M 2 (we refer to this 
approximation as the "OPE-limit"), can be found in ref. [14]. Since our techniques work 
beyond the OPE-limit, we expect them to allow for an improvement of such computations. 
A few other cosmological applications are suggested in the Conclusions. 

Previous work 

The very first spectral function computed at next-to-leading order (NLO) in the domain 
Lxj ~ 7rT, p = was that related to the electromagnetic current associated with massless 
quarks [15]-[17]. These results have subsequently been generalized to the case of a finite 
(heavy) quark mass [18], and another spectral function related to this case, correlating Lorentz 
forces felt by the nearly static heavy quarks, was derived at NLO with similar methods [19]. A 
general analysis of many spectral functions in the domain uj 3> irT, which can be systematized 
through Operator Product Expansion (OPE) techniques, was carried out in ref. [20]. The 
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current study aims to generalize the approach presented in refs. [18, 19]. Note that there is 
a huge body of literature concerning the domain oj < gT, p = (here g = ^4ira s ), in which 
perturbation theory breaks down; upon this we touch only briefly in sees. 5.2, 5.3. 

Outline of this paper 

The observables that we are specifically concerned with are the 2-point correlators of the 
scalar ("FF") and pseudoscalar ("FF") densities in pure SU(iV c ) Yang- Mills theory. These 
quantities have recently been studied particularly by Meyer, who has addressed both the 
Minkowskian and Euclidean domains [21]-[24]. Noteworthy is also a general discussion con- 
cerning the shape of the scalar channel spectral function at small frequencies [25] as well 
as the mentioned OPE representation of its ultraviolet asymptotics [20]. In addition, sum 
rules have been analyzed [26], improving on earlier findings [27, 28], and further insight has 
been sought from AdS/CFT inspired models [23, 29, 30, 31]. The corresponding transport 
coefficients (this yields the bulk viscosity for FF, cf. e.g. ref. [32], and the rate of anomalous 
chirality changing transitions for FF, cf. ref. [33]), together with possible applications, have 
been the subject of much further work; we restrict here to mentioning determinations of the 
transport coefficients in the weak-coupling limit [34, 35]. 

The current work is a continuation of two recent papers which also addressed the same 
general problems. In ref. [36], the scalar and pseudoscalar correlators were studied within 
the OPE framework, refining the previous determination of certain Wilson coefficients and 
also correcting results concerning the pseudoscalar channel [23]. In ref. [37], which went 
beyond the OPE-limit, the full |x|-dependences of the time-averaged Euclidean correlators 
were determined, and compared with related lattice measurements carried out in ref. [23]. 
One specific purpose of the present paper is to "complement" ref. [37], by generalizing the 
results beyond the OPE-limit also in the Minkowskian frequency domain. 

As already mentioned our current analysis borrows techniques from refs. [18, 19] but goes 
beyond these works in its generality. We hope that we can thereby pave the way for a num- 
ber of further generalizations, such as the inclusion of the Lorentz non-invariant structures 
appearing e.g. in the "shear" channel of the energy-momentum tensor; of a non-zero spatial 
momentum, relevant for relativistic particle production rates in cosmology; and of effects 
from a non-trivial mass spectrum, relevant e.g. in chiral perturbation theory applications. 

The paper is organized as follows. The observables are defined in sec. 2, the general method 
for their determination is discussed in sec. 3, and the basic results are presented in sec. 4. 
The results are analyzed in some detail in sec. 5, which also contains numerical evaluations. 
We conclude and offer an outlook in sec. 6. Three appendices, A, B and C, contain a detailed 
investigation of one of the sum- integrals, concise results for the others, and an account of a 
resummation relevant for g 2 T/ir <C co <C vrT, respectively. 
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2. Specific setup 



Our basic notation follows ref. [36], so we discuss the specific setup only briefly. We consider 
pure SU(-/V C ) Yang-Mills theory (with N c = 3 in numerical estimates), dimensionally regular- 
ized by analytically continuing to D = 4 — 2e space-time dimensions (the Euclidean action 
reads Se = Jo dr fd 3 ~ 2e XjF^ u F® u ). The operators considered are 

where g\ is the bare gauge coupling squared; these operators require no renormalization at 
the order of our computation. We normally leave the coefficients c 9 , c x unspecified, but note 
that often the values c d pa — ^ — , c x = are chosen, where b = and b\ = ■ 

The bare gauge coupling squared can be expanded in terms of the renormalized one, g 2 , as 



2 2 2e 

9b = 9 V 



(2.2) 



where \i denotes a scale parameter. The MS scheme renormalization scale is denoted by 
p? = 47r / u 2 e~ 7E . In order to avoid unnecessary clutter all appearances of /U 2e , which play no 
role in the final renormalized results, will be suppressed. 

The observables considered are related to the 2-point correlators of 9 and x- We define the 
quantities (X = (r, x); P = (p n , p); p n = 2irTn, n G 7L) 

G e (X) = ( 6{X) 6(0) ) T , G X (X) = ( X (X) x (0) ) T , (2.3) 

as well as the corresponding Fourier transforms 

Go(P) ^ [ e- ip - x G e (X) , G X (P) ^ [ e~ ip - x G x (X) , (2.4) 

where J x = J^dr J x . No hats appear above the operators because Euclidean correlators can 
be evaluated with standard path integral techniques. It can be shown (cf. e.g. refs. [1, 2]) 
that the corresponding spectral functions, defined as in eq. (1.1) with O — > 9,x, are then 
obtained (at vanishing spatial momentum) from 

p(u) = Im\ G(P)] . (2.5) 

It is these quantities (denoted by p , p x ) that we are mostly interested in. 



3. General method 
3.1. Wick contractions 

Having defined the correlators, the first step in their determination is to carry out the Wick 
contractions (for our specific example the corresponding graphs are shown in fig. 1), and to 
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(ii) (iii) (iv) 

O <x> <X> 

(v) (vi) (vii) 

Figure 1: The graphs contributing to the 2-point correlators up to 2-loop order. The wiggly lines 
denote gluons; the small dots the operators 8 or \ ( CI - e q- (2-1)); and the grey blob the 1-loop gauge 
field self-energy. "Disconnected" contractions only arise at higher orders. 

reduce the functions Gq{P), G x {P) to sums over various types of integrals. We strongly 
advocate taking this step in Euclidean signature, whereby the Feynman rules are elementary 
(no excessive i's or doublings of degrees of freedom appear), and by using bare parameters, 
whereby counterterm vertices are avoided. An additional strength of the imaginary-time 
formalism is that in terms of physical processes, all "real" and "virtual" corrections related 
to each other are automatically captured by a single Matsubara sum-integral. Before carrying 
out the contractions, it is helpful to make use of translational invariance in order to write 

G e (P) = [ e- ip - x (9(X)8(0)) T = [ e^ p < x ~ Y \e{X)9{Y)) T 
Jx Jx 

= j-(J x e- ip - x 9(X) j Y e^ Y 9(Y)) T , (3.1) 

so that the full analysis can be carried out in momentum space (the factor J Y = /3V = 
j38 pn=0 (2ir) d 5( d Xp = 0), with d = D - 1, cancels out because of momentum conservation). 

3.2. Scalarization 

In the next step, the goal is to "scalarize" the sum-integrals, i.e. to turn them into ones that 
also appear in scalar field theories. This can be achieved by contracting Lorentz indices (using 
= D) and carrying out Dirac traces (with rules following from {7^, 7^} = 25^ u l 4 x 4 ; in the 
case of chiral fermions the usual issues with 75 need to be faced and we have nothing to add on 
this topic). Subsequently, substitutions of sum- integration variables, completions of squares 
(e.g. Q ■ R = 7} [Q 2 + R 2 — (Q — R) 2 ]), and identities following from partial integrations with 
respect to spatial momenta, can be employed in order to remove as many scalar products as 
possible from the numerators. The goal is to express the result in terms of a minimal number 
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of independent "master" sum- integrals. (If the set is minimal indeed, then all dependence on 
the gauge fixing parameter must have cancelled after this step.) 

For the specific correlators Gq(P) and G X (P), the reduction to master sum-integrals was 
carried out already in ref. [36]. Re-expressing a sum-integral "X;" of ref. [36] in terms of "Xj," 
(cf. eq. (A.14) there) the results from eqs. (3.1), (3.2) of ref. [36] become (d A = N 2 - 1) 

G e {P) 



34- 13D 
+ 5 X i~ 



-(D-2)X a + (D-4)X b 
<D-A) 2 



-l R + (D-2) 



+ {D- 2) 2 
-X e + 3X h + 2X i ,-X j 



(3.2) 



G X (P) 



-l6d A c 2 (D - 3) 



9i(D-2) 



+ glN c \2{D-2){D-A)l h + {D-2) 



2D 2 - IID + 30 



X f - 2{D - 4)X„ + (D — 2) 



-X c + 3X h + 2X i ,-X j 



. (3.3) 



Here, for brevity, structures containing $ Q 1, which vanishes exactly in dimensional regular 
ization, have been omitted = TY^ qn J* q )- The master sum-integrals are defined as 

p2 

,Q 2 ' 

p4 



UP) 
MP) 
UP) 
UP) 
UP) 
UP) 
UP) 
UP) 



Q 2 (Q-P) 2 



Q,R 



Q,R 



Q,R 



Q,R 



Q,R 



Q,R 



Q 2 R 2 ' 

p2 


Q 2 R 2 (R- 


P) 2 


p2 




Q 2 R 4 ' 




p4 




Q 2 R 4 (R- 

p2 


P) 2 



Q 2 R 2 (Q - R) 2 ' 



Q 2 (Q-R) 2 (R-P) 2 



(3.4) 
(3.5) 
(3.6) 
(3.7) 
(3.8) 
(3.9) 
(3.10) 
(3.11) 
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Xg(P) ~ T QR Q 2 (Q-P) 2 R 2 (R-P) 2 ' (3 ' 12) 



lh{P) T Q 2 R 2 (Q-R) 2 {R-P) 2 ' (3 ' 13) 

Q,R 

MP) 55 $ n &R*(Q-R)*(R-P)*> (3 ' 14) 

Xj(P) = ^ Q 2 i? 2 (Q - R) 2 (Q - P) 2 {R - P) 2 " (3 ' 15) 
3.3. Matsubara sums and discontinuities 

The master sum-integrals denned in eqs. (3.4)-(3.15) contain a 1- or 2-fold Matsubara sum. 
Both sums can be carried out explicitly. In the literature various recipes can be found for 
this, most notably "cutting rules"; we have carried out the sums with the method described 
in some detail in appendix A.l of ref. [18] and illustrated with simple examples also in 
eqs. (C.13)-(C.15) and (C.29)-(C.30) below. Partial fractioning the result in the variable p n , 
there may be polynomial terms, as well as fractions of the type 

(3.16) 



where = ±1 and 6 {-E q , -Eq-p, E r , E r _ p , -E q _ r }, E^ = |q|. According to eq. (2.5), the 
corresponding spectral function follows from 



Im 



nqr-U =- =-7r6(u + J2<7kEk). (3-17) 



The spectral functions corresponding to all the master sum-integrals after this step (omitting 
terms oc uj n 5(ui) after setting p — > 0, cf. sec. 5.6) are listed in appendices A and B. 



3.4. Spatial integrals 

The "hard work" of the problem, where little automatization seems possible, is to carry 
out the remaining spatial integrals, constrained by the Dirac-J's. There are two classes of 
integrals, referred to in the literature as "virtual" and "real" corrections; we often refer to 
them as "factorized" and "phase space" integrals, respectively. Within the former class, 
some integrals are ultraviolet divergent, necessitating a careful handling by e.g. dimensional 
regularization. All classes of integrals are described in detail in appendix A for the specific 
spectral function denoted by p x , corresponding to the sum-integral defined in eq. (3.15), and 
more concise results are collected in appendix B for all the other cases. 
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3.5. Renormalization 



After all master spectral functions are known, these are inserted into (the imaginary parts 
of) eqs. (3.2), (3.3). Simultaneously, D = 4 — 2e is inserted into the coefficients, and the bare 
gauge coupling is re-expanded in terms of the renormalized one according to eq. (2.2). If the 
computation has been carried out correctly, all 1/e-divergences must cancel at this stage, and 
we can subsequently set e — > 0. 

Let us illustrate the procedure with the example at hand. A number of the sum-integrals 
defined have no imaginary part X a , X c , X c ), and therefore lead to a vanishing spectral 
function. Among the non-zero spectral functions, only three have an 1/e-divergence (p x , 
p x , p x ). Taking these facts into account, eqs. (2.5), (3.2), (3.3) can be converted into 

e) {gi p Jh M + 2g% N c [3^ (oj) + 2p x , (oj)] } 

- 2glN c [2p Xd (u) + 3 tef M + Pi .(oj)] + 0(g e B e,g 8 B ) , (3.18) 

2e)(l-e){gi Pjb (u) + 2glN c [3 PXh (u) + 2p Xi Xu>)]} 

- 2gtN c [2p Xd (u) + 3/> Xf M + Px .(oj) - 2ep Xg (oj)] + 0(g 6 B e,g 8 B ) . 

(3.19) 

The structures within the curly brackets on the first rows turn out finite after the insertion 
of the bare gauge coupling from eq. (2.2) and the re-expansion of the result in terms of the 
renormalized coupling; therefore the prefactors 1 — e, 1 — 2e can actually be set to unity. 

3.6. Final expression and its limits 

After renormalization and setting e — > as outlined in sec. 3.5, we obtain our final results 
for the spectral functions. Furthermore, if they are only needed in the OPE-regime oj S> ttT, 
then all master spectral functions can be determined in a closed form (cf. ref. [36] for general 
methods as well as appendices A and B for specific results for the structures in eqs. (3.4)- 
(3.15)). Moreover, the general outcome can be given an interpretation in terms of thermal 
contributions to "condensates" , cf. ref. [20] , which guarantees the cancellation of terms oc T 2 , 
if the lowest-dimensional gauge-invariant condensate has dimensionality four, as is the case 
in our study. (If scalar fields are present, then there can be a contribution proportional to 
the condensate ~ (<fy<f>), implying that terms oc T 2 need not cancel; see e.g. ref. [14]). 

The situation is considerably more subtle in the infrared regime, oj ttT. In fact, at 
small frequencies, oj ~ gT, the NLO corrections become as large as the leading-order (LO) 
terms, indicating a breakdown of the perturbative series; this will be discussed with our 
specific examples in sec. 5.2. In this situation perturbation theory needs to be resummed 
through effective field theory techniques. For oj ~ gT the relevant framework is that of the 
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Hard Thermal Loop effective theory (ref. [38] and references therein), whereas for smaller 
frequencies still, needed e.g. for transport coefficients, a further re-organization is required 
(ref. [34] and references therein) . The different infrared frequency scales playing a role in the 
scalar channel have been elaborated upon in ref. [25]. 



Results 



Inserting the spectral functions pj , p x , p x from eqs. (B.6), (B.44), (B.56), respectively, as 
well as the bare gauge coupling from eq. (2.2), all divergences are seen to cancel in eqs. (3.18), 
(3.19), as already advertised. Making also use of the finite functions p x , p x , p x from 

d f j 

eqs. (B.19), (B.27), (A. 57), respectively, as well as the divergent part of p x from eq. (B.33), 
the final results can be expressed as 



Mac 2 , 



(l + 2n«)<g + 



-P x \ u ) 
IQcLac 2 , 



(47T) 



TTOJ 



(l + 2n f )i«7 4 + 



(4tt) 2 [ 3 

9 6 N< 
(4vr)' 



22 ln£ 



73 

+ + 



22 pf- 
— In — ^ 

3 iu 2 



97 



\4> T (u) 



+ 0(g 8 ) 
+ 0(g 8 ) 



(4.1) 
(4.2) 



(4vr) 2 

where we have introduced 

n E = n B (E) , (4.3) 

with n B (E) = \j[eP E — 1) denoting the Bose distribution. To this order the only differ- 
ence between the two channels is in the constant next to the logarithm. The logarithms 
and the constants next to them agree with those found within the OPE-regime in ref. [36]. 
(In principle they were computed already in ref. [39], however in a different regularization 
scheme.) The "non-trivial" T-dependence resides inside the function <f> T , which is given by 
(we have substituted q = ua/2, r = ujt/2 in the expressions of appendix A and B, and 
defined h x = n<^) 



0tM 



-L 



d<7 h n 



+ 



1 1 

a a — 
1 

a a 



-2 + a-a z 



+ 2 + a + a 2 



ln(l - a) 
ln(l + (j) 



11 

+ 12 



1 



+ 



/: 



da h a 



CT + 1 

2 



+ 



1 



1 

a (j — 1 

"1 1 

(7 (7 + 1 



11 



5(7 

+ T 

2 



+ 



+ 



1 



(7 



3 (7 

T + 4 + ^ 



+ 2 + (7 + (7^ 
2" 



ln(l - a) 
ln(l + (r) 



In(a) 
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+ 



da h c 



+ 
2 
+ 
+ 
+ 



11 1 1 (7 

2a 

12(7 + 1 3 3 



1 



1 



a a — 1 
1 1 



2 + a - a 1 



+ 2 + + CT 1 



a a + 1 

1 1 5 

- - + 2- ct + <t 2 

a — 1 <t 



In (a - 1) 
ln(l + (j) 
ln(a) 



11 

12 



+ 



o jo 



1 



1 



a + 1 a 



23 13a lla 2 
_ _ _ _ 
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2 l CT 2I ni_ CT n CT+T (l +ni_ 
dr m 



X 



X 



1 5 - 4r + 2r 2 
ar 4a 



5 - 4cr + 2a 2 3 
4r" + 2 



+ 



00 r<J—l 

da / 

Jo 



+ 



X 

da / dr 
Jo 

X 



n CT _i(l + n a _ T )(h a - h T+1 ) 

CIT ; — X 

n T n_i 

1 5 + 4r + 2r 2 5 - Aa + 2a 2 3 
ar 4cj 4r 2 



(1 +n a+l )n a+T h T+1 



n 2 

1 5 + 4r + 2r 2 5 + Aa + 2a 2 3 

— + 'A + A + O 

ar 4a 4r 2 



(4.4) 



It is conceivable that some of the structures in eq. (4.4) could be written in a simpler 
form; given that such rewritings do not allow to reduce the dimensionality of the integration, 
however, we prefer to display the result in the current form, suitable for numerical evaluation. 
On this point we should mention that eq. (4.4) is defined in the sense of principal value inte- 
gration; this implies that the 2nd terms from the 2nd and 3rd structures could be combined 
into 

/ | dfT ^{-^T ln ( 1 - CT ) } + / da^j-^-ln^-l) I 



(4.5) 



or, alternatively, that the range (1, 00) in the 3rd structure could be reflected to (1, 0) through 
a — > 1/a and then combined with the 1st and 2nd structures. In any case the integral is 
rapidly convergent; the result of eq. (4.4) is illustrated in fig. 2. 
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Figure 2: The function (j> T from eq. (4.4), multiplied with (^) 3 , as a function of -tW. Also shown 
are the ultraviolet ("OPE") limit from eq. (5.1) as well as the infrared ("IR") limit from eq. (5.10). 

5. Analysis of the result 
5.1. Ultraviolet limit (to > 7rT) 

As has already been mentioned the general results simplify considerably in the OPE-limit 
to S> vrT and can be expressed in a closed form. Indeed, according to refs. [20, 36], the 
asymptotic to 3> vrT behaviour of the function <j) T from eq. (4.4) should read 



Let as stress again that, as required by OPE [20], all contributions of 0(tt 2 T 2 /to 2 ) from 
individual graphs need to cancel in the sum, because no gauge-invariant condensate of di- 
mensionality GeV 2 is available. It is easy to verify numerically that the form of eq. (5.1) is 
indeed approached, cf. fig. 2. The power like asymptotic behaviour comes solely from the first 
integral in eq. (4.4), the other terms falling off exponentially at w > 7rT. 

5.2. Infrared limit (^T/vr < to < ttT) 

Another interesting limit is the infrared one, to <C ttT, which shows a much richer structure 
than the ultraviolet limit. Denoting the six separate parts of eq. (4.4) by Jj, . . . , I e , respec- 
tively, with the combination appearing in eq. (4.5) taken apart and denoted by I^, some 




(5.1) 
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work leads to the small-w expansions 

Ji = {-2Li 2 (|) -21n3arcoth7-^^ -ln 2 3 + 61n2+4 -|}^ , (5-2) 

1 2 = |Li 2 (i) +21n3arcoth7+^^ +31n 2 2-81n2 + ^ -3}^ +0(w) , (5.3) 

1 3 = - 44C 3 ( ^ r3 + { -8 ln(ft ) - 96 ]nA + 8 ln(47r) + 7} ^ 

+ {-f ln(| ) - ln 2 2 + M^ + iiV + 37j^ + 0(]nw) _ (54) 

J 2 d ] v = f2Li 2 (i) -4Li 2 (§) -21n 2 3 + 41n 2 2- ^j- +0(lnw) , (5.5) 



T 



{¥-T+n}-+OM, ( 5 - 6 ) 

2tt 2 T 2 

J 6 = |ln(|) + 121nA-ln(4vr)-|}^-^ 

+ {ln 2 (^)-(21n2 + i)ln(^) + ci|- +0(hiw) , (5.7) 

L > id 

h = 11 |P+{l^) + 121n^-ln(4 7 r) + ^}^ 

+ {-ln 2 (^) + (21n2 + f )ln(^)+c 2 }^ +0(lnu) , (5.8) 

where A stands for the Glaisher constant, In A = — £(— 1) — £'(— 1), and we have defined the 
numerical coefficients 

ci » -0.17449 , c 2 » -9.32085 . (5.9) 

Summing together, terms of 0(T 3 /u; 3 ) as well as all the logarithms cancel, and we are left 
with the expansion 

(f) T (oj) "%" T ^ + 3.31612 ^Wln^). (5.10) 
6L0 UJ V 7rT/ 

As can be seen from fig. 2, this does agree with a numerical evaluation of 4> T for oj <C vrT. 

Now, even though we can derive a representation of the function 4> T for oj <C 7rT, this does 
not mean that we know the spectral function in this regime. The reason is that, as can be 
seen from eqs. (4.1), (4.2) and (5.10), the thermal NLO corrections overtake the LO terms 
for u <C gT. Therefore the naive perturbative expansion breaks down in this regime. The 
spectral function can be determined only if we find a resummation which re-organizes the 
perturbative series in a way that a breakdown is avoided. 

It is strongly believed that the way to resum the perturbative series around oj ~ gT goes 
through the use of the Hard Thermal Loop (HTL) effective theory (ref. [38] and references 
therein). Although caveats may be difficult to exclude [40], we demonstrate in the following 
that this framework does indeed reproduce the leading divergence of eq. (5.10) and, through 
a matching computation, should allow us to postpone the breakdown tow~ g 2 T '/ir. 
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Schematically, matching can be represented as 



.QCD 
r resummed 



QCD _ HTL , HTL 

r resummed r resummed 1 i resummed 



QCD _ HTL , HTL 

r naive r naive 1 r resummed 



(5.11) 



Here we have subtracted and added the HTL part, and subsequently noted that if the correct 
effective theory is used, then the difference between the full and the effective theory compu- 
tations is infrared safe, so that no resummation is needed for the difference. Since in sec. 4 a 
"naive" QCD computation was reported, we then need two different HTL computations, one 
naive and one resummed, in order to obtain the correctly resummed version for QCD. 

Some details concerning the HTL computations are given in appendix C. As far as the 
naive version goes we note that, replacing the Bose distributions through their "classical" 
limits, i.e. n q — > T/q, 1 + n u _ q — > T/(ui — q), etc., which is equivalent to the physics of HTL 
resummation [41, 42], the integrals in eq. (C.25) of appendix C.2 are elementary and we get 



4d A ci 



16dACy 



-Kg 



4 (l + 2nt 



(47T) 



V + aA^}+0(<? 8 ). 



(5.12) 



Here the Debye mass parameter is defined by 



g 2 N c T? 



(5.13) 



The first term of eq. (5.12) matches the leading-order QCD result in eqs. (4.1), (4.2), whereas 
the second term exactly corresponds to the leading term of eq. (5.10). So, in the difference 



QCD 
A^naivc 



Pnlivc appearing in eq. (5.11) the dominant IR divergence drops out. 



As far as the resummed version of the HTL computation goes, it is technically similar to 
the classic dilepton analysis in ref. [43]. The main results are given in eqs. (C.49), (C.50) of 
appendix C.3. Let us define 



4d A ci 



16gU<4 



resummed 



resummed 



7rg 4 (l + 2n f ) 

(4vr) 2 
7rg 4 (l + 2n f ) 

(4vr) 2 



oj a + u) 2 m\ + m 4 



UJ 



</>r L M} 



HTL 

X 



(5.14) 
(5.15) 



so that the separated terms cancel against the naive version in the difference appearing in 
eq. (5.11). The result is shown in fig. 3. (The first "cusp" corresponds physically e.g. to a 
2o2 process in which one on-shell "plasmon" gets generated; the second to the 2-plasmon 
threshold.) For u < m E , </>#^ L modifies the IR behaviour of the full expression in a qualitative 
way (cf. fig. 4 below). For co » m E , in contrast, (j^ L is a constant of 0(g A m E ) ~ 0(g 8 T 4 ) 
which is vastly overshadowed by the unresummed behaviour ~ g 4 co 4 ; therefore in this regime, 
which is the most important one for us, HTL resummation plays no role. 
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Figure 3: The functions <pg T, L from eqs. (5.14), (5.15), multiplied with ^E. ( so th a t the same number 
of overall ui's has been factored out as in fig. 2). Also shown are partial contributions from various 
sub-processes ("-naive" refers to the terms w 4 + Lo 2 m\ taken apart as in eqs. (5.14), (5.15); in the left 
panel, thick lines correspond to "transverse" and thin to "electric" modes, respectively). 



5.3. Extreme infrared limit (lo < g 2 T/n) 

The HTL resummation discussed in sec. 5.2 is supposed to re-organize the perturbative series 
in a way that it can formally be defined even in the range u ~ gT. If, however, the frequency 
is decreased further, down to to < g 2 T/-ir, then a further re-organization is needed in order to 
obtain the correct leading-order result [25]. The "transport peak", whose height yields the 
transport coefficient, appears within this regime; in our computation part of the information 
concerning it is hidden in the so-far omitted structures ~ co n 5(uj), to be discussed in sec. 5.6. 



5.4. Sum rules 

Apart from the UV and IR limits, an interesting crosscheck on the full result is offered by 
"sum rules", which concern integrals over the spectral function with a certain weight. These 
integrals are related to "thermodynamic" Euclidean observables, which can be evaluated 
independently. Let us elaborate on the information that can be obtained this way. 

Because of issues of convergence, some terms need to be subtracted before sum rules can 
be defined. If we formally subtract the zero-temperature parts from both sides, denoting 
the results by "A" (in practice such a subtraction is not without problems, as discussed in 
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sec. 5.5 below) then, to the order of our computation, the sum rules boil down to [36] 



ff 6 iV c r 4 _ AG e (0) 
Ad A c 2 9 

AG x (0) 



18 



-16d A c x 



Jo ™ 
Jo iru 



4d A c 2 
-Ap x (a; 
16dAc| 



+ 0(5 8 ) 



(5.16) 
(5.17) 



Here, the thermal parts read (from eqs. (4.1), (4.2)) 



Ap e {u) 


4 

7T0J 


4d A c 2 e 


(4vr) 2 


-Ap x (w) 
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(4tt) 2 
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2rii 
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(4tt) 5 
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+ (l + 2n«) 



4 5 b ^c 

<? 4 + 



(47T) 



22, u 2 97 
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+ (l + 2n<-)- 



8g*N c <t> T (u) 
(4vr) 2 

fN c cf> T {uj) 



(4vr) 2 



(5.18) 



(5.19) 



Now an immediate "paradox", raised in ref. [25], is that even though eqs. (5.18), (5.19) 
have a common positive-definite thermal part at 0(g 4 ), the integral over this is supposed to 
vanish, since the left-hand sides of eqs. (5.16), (5.17) are of 0(g 6 ) or higher. A resolution 
to this paradox in the renormalized case, in which the coupling constant runs, was given in 
ref. [20]. Beautiful as the argument is, underlining the strength of sum rules in that they 
allow us to anticipate features of a higher-order computation without carrying it out, the 
method is cumbersome in our case when we want to use the sum rules as an exact crosscheck 
of an order-by-order computation. In the following we therefore "freeze" the gauge coupling, 
so that different orders do not mix. The price to pay is that then the issue of "contact 
terms", discussed e.g. in ref. [26], needs to be addressed. On the other hand, the contact 
terms necessarily play a role in the "shear" channel [20, 44], where similar structures appear 
but without g A and running, so the discussion may be useful as an analogue. 

Contact terms arise when there is a part in the Euclidean correlators, AGe(P) and AG X (P), 
which does not vanish as pl^ — > oo. This part is "lost" to the spectral function, and needs to be 
"added" to the right-hand side of the sum rule, before a comparison with the thermodynamic 
left-hand side, in which the information is included, can be carried out. (A more precise 
account of the logic can be found, e.g., in ref. [26].) 

From eqs. (4.1) and (4.2) of ref. [36], setting P = (p n , 0) and inserting J^qn B (q) = 7r 2 T 4 /30, 
Jq n B (q)/q = T 2 / 12, the ultraviolet limits read 



AGe(P) Pn »nT 47ryT 4 



15 



M A c 2 e 

AG X (P) Pn> >, T i^g^T 4 



-16d A c x 



15 



1 + 
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g 2 N c 

(4vr) 2 
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(4vr) 2 



22, fi 2 203 



22, a 2 347 
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' >N cT 4 „ 
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Pr, 



18 
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(5.20) 



(5.21) 
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where we have kept separate the terms coupling to two different operators in the OPE limit. 

At 0{g i ), we now directly observe that if the contact contributions from eqs. (5.20), (5.21) 
are added to the integrals over the leading-order terms in eqs. (5.18), (5.19), which evaluate 
to 

2 9 / — 7T^2 2n f = "T" / dxx n * = k > 5 - 22 

J 1TUJ (AttY - ; IT 2 Jq 15 

then the desired cancellation duly takes place. 

As far as 0{g & ) goes, we observe that the constants next to the logarithms in eqs. (5.20), 
(5.21) can be written as^ = f- ^f,^f = f- ^f. By making use of eq. (5.22), the 
parts containing ^ and ^ are seen to exactly cancel against corresponding terms from the 
integrals over the spectral functions in eqs. (5.18), (5.19). Therefore, the non-trivial sum rule 
has a common — ^ in both channels. We also observe that the difference in the last terms 
of eqs. (5.20) and (5.21) precisely explains the difference of the left-hand sides of eqs. (5.16), 
(5.17). So, we are left with one non-trivial sum rule to check; we choose eq. (5.17). 

It requires some care to formulate this sum rule. The reason is the asymptotics of eq. (5.1) 
which, for a frozen coupling, implies that the integral in eq. (5.17) diverges. This is reflected 
by the appearance of the logarithms involving p 2 on the right-hand side of eq. (5.21) (the 
imaginary part of this logarithm gives the asymptotics of eq. (5.1)). Let us regulate the 
divergence by restricting u to \oj\ < A. Noting that 

(1 \ f A du; 
hw = 7T , 2/ — 7r = lnA 2 , (5.23) 

we deduce that in the contact term we can replace In p 2 — > In A 2 in the presence of the cutoff. 
With this logic, the coefficients of g 6 N c in eqs. (5.17), (5.19) and (5.21) imply the sum rule 
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It)- <5 - 24) 

where — ^ = — ^ — . We have verified numerically (with a relative error smaller than 
3^5 ) that this sum rule is satisfied by our function <p T . 

Sum rules could also be used in connection with the HTL-results in eqs. (5.14), (5.15). 
As far as the thermodynamic left-hand sides are concerned, eq. (C.7) yields a contribution 
3^(^)2 g 7 T A to eq. (5.16), which equals the next-to-leading order correction to T 5 ^y(^|£), 
whereas eq. (C.8) implies that the left-hand side of eq. (5.17) vanishes even in the presence of 
HTL-resummation. However, there are again contact terms which need to be added in order 
to satisfy the sum rules. Given that HTL-resummation is of secondary importance for our 
study, we refrain from a more detailed discussion here, noting only that a cancellation in the 
X-channel is also visually suggested by fig. 3(right). 
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5.5. Numerical evaluation 

After the various crosschecks carried out in sees. 5.1, 5.2 and 5.4, we now move on to the 
numerical evaluation of our results. The goal is to obtain curves which can in principle be 
used in the applications outlined in the introduction. 

In order to evaluate our result (eqs. (4.1) + (4.2) + (5.11)) numerically, we need to assign 
a value to the running coupling g 2 . It is only in the regime u 3> ttT, where 4> T <C 1, 
that two subsequent orders with the same functional form are at our disposal, so that a 
scale optimization is possible; we then define Pq P ^ from the "fastest apparent convergence" 
criterion based on eqs. (4.1), (4.2): 

ln(/2° pt(w) )^lnH-g, ln( A fH)^lnH-g. (5.25) 

In the infrared regime u <C ttT we choose g 2 from the "EQCD" setup (cf. ref. [45] and 
references therein) which at 1-loop order amounts to: 

ln(/2° P x t(T) )^ln(4^T)- 7E -l. (5.26) 

For a given oj the larger among these scales is chosen; the switch happens at oj ~ 11.3 irT for 
p g , and at oj ~ 19.5 ttT for p x . We note that the extremely "late" transition to the vacuum 
scaling is related to the well-known slow convergence of the perturbative expansion for the 
Wilson coefficients appearing in the OPE regime. Indeed the fact that in the full computation 
there is another scale present, ttT, which acts as a sort of an infrared cutoff, allows us to 
obtain a numerical evaluation for frequencies much smaller than in the OPE regime [36]. 
There is a certain price to pay, however, which is that it is no longer possible for us to split 
the result into a "zero-temperature" and a "finite-temperature" part, because they have been 
mixed by the scale choice. (Put another way, the two parts separately are quite sensitive to 
the confinement scale, but their sum, which we evaluate, is not.) 

In order to get an error band for the uncertainty related to the choice of g 2 , we use the 
2-loop running of g 2 , with the renormalization scale varied in the range (0.5. .. 2.0) x p opt . 
The scale parameter, defined as A^g = lim^oo p [bog 2 ] bl ^ 2b ° exp[— — Uj ], is re-expressed 
through the critical temperature of the deconfinement transition of pure SU(3) gauge theory; 
we take this to be T c = 1.25A^g-, which could also be viewed as our definition of "T c ". 

The outcome of a numerical evaluation of eqs. (4.1), (4.2), (5.11) is plotted in fig. 4, in 
units of uj 2 (ttT) 2 . We observe that the dependence on the scale choice practically disappears 
from the NLO result in the ultraviolet domain, oj S> vrT. As far as the infrared domain is 
concerned, we observe an enhancement due to HTL-resummation in the ^-channel akin to 
that seen in the vector channel [43] , whereas in the x-channel the HTL-resummation removes 
much of the spectral weight from small frequencies. 
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Figure 4: A numerical evaluation of eqs. (4.1), (4.2), (5.11), in units of w 2 (7rT) 2 , for T = 3.75A|3g 
corresponding to T = 3T C . The gauge coupling has been fixed as explained around eqs. (5.25), (5.26), 
and the gray band reflects the corresponding uncertainty. (In the 0(g 4 ) result the "optimal" scale is 
fixed to the thermal value of eq. (5.26), i.e. does not change with the frequency.) 



In ref. [24], a significant "cancellation" in the thermal part of the spectral function was 
advocated for, with a positive domain at small frequencies and a negative domain at large 
frequencies. For reasons mentioned above, we do not believe our numerical evaluations to be 
quantitatively reliable for the zero-temperature and finite-temperature parts separately, and 
therefore abstain from a direct comparison. Nevertheless, we note that on the qualitative 
level the structure of the finite-temperature part as discussed in ref. [24] (cf. fig. 5 there) 
agrees perfectly with that of the function <j) T as shown in fig. 2. 

5.6. Imaginary-time correlators 

Although the ultimate goal of a lattice study could be the determination of a transport 
coefficient, which is by definition a Minkowskian object and obtained from the zero- frequency 
intercept of p(u,0)/uj, any practical measurement is carried out in Euclidean signature. It is 
relatively easy to go from a known Minkowskian spectral function to a Euclidean correlator, 
but very difficult to take the inverse step; therefore, as a first comparison between lattice 
data and an analytic computation, it may be helpful to extract the Euclidean correlator from 
the latter. The principal tool for this is the well-known formula 

roo j COsh ( %• — T ) U 

G{r)= -p(u,0) Y ' . (5.27) 

Jo 71" smb. 
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As it turns out, there are some problems hidden in eq. (5.27). One is the issue of "contact 
terms", discussed in sec. 5.4; as the name says, one may move them to the Euclidean side 
of the relation and interpret as a zero-distance term oc <5(r), giving a finite contribution to 
Jq 'dr G{t). Here we stay at non-zero r, and this issue does not arise. However, there is 
another issue, namely the possible existence of terms of the type ~ u)d(u>) in the spectral 
function which, according to eq. (5.27) (replacing f£°duj — > \ dco) can give a constant 
contribution to the Euclidean correlator. This can be important particularly around the 
middle of the Euclidean time interval where the absolute value of the correlator is smallest 
(cf. e.g. ref. [46] and fig. 5 below). 

Before proceeding with the specific discussion, let us note that in general the two issues 
mentioned are somewhat "complementary". Operators with a higher dimension are more 
ultraviolet sensitive; this means that, like in many of the sum-integrals in eqs. (3.4)-(3.15), 
explicit powers of P 2 appear in the numerator and therefore contact contributions, originating 
from a non-trivial limit of the Euclidean correlators at p 2 t 3> (ttT) 2 , play a role; cf. sec. 5.4. If 
the operator has a lower dimension, with less or no powers of P 2 , then the ultraviolet regime 
is safer; on the other hand, the issue with zero-mode contributions is more likely to arise. 
(Only in exceptional cases are both problems avoided [47, 19].) 

Let us recall a simple illustration of the issue. Consider the structure defined in eq. (3.5), 
but without the P A in the numerator, 

Carrying out the Matsubara sum and taking the imaginary part leads to (we keep p / for 
the moment; E qp = |q — p|) 



f—{ 

Jq AqE qp { 



5(u - q- E qp ) - 5(uj + q + E qp ) (1 + n q + n qp . 



+ 



8(u> - q + E qp ) - 5{uj + q- E qp ) (n, 



"qp 



(5.29) 



(This is like eq. (B.4) but now with the second term added.) For oj — > the first row gives 
no contribution because of vanishing phase space, whereas on the second row we can write 



5(u> - q + E qp ) - 5(u + q - E qp ) (n qp - n q ) 



S(uj - q + E qp ){ 



n, 



q—u 



n, 



,) - S(u + q - E qp ) {n q+u) - n q ) 



— —2uS(uj)n' = 2/3uj5(uj) n q (l + n q ) . 

Inserting this into eq. (5.27) (with / °°dw — > \ J^dco) yields 

/oo r ^ 2 f l 



(5.30) 



(5.31) 
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Figure 5: A numerical evaluation of Ge (left) and G x (right), in units of T 5 , for T = 3.75Aj^g 
corresponding to T = 3T C . The gauge coupling has been fixed as explained in fig. 4, and the gray 
band reflects the corresponding uncertainty. 



which would yield a constant (T-independent) contribution of the type advertised (in this 
case it is even infrared divergent). 

In our situation, however, eq. (5.28) contains an additional P 4 or, after analytic contin- 
uation, uj 4 . This kills the contribution from u8{uS) at 0(g 4 ). At higher orders, the loop 
expansion may introduce additional factors of g^n^ ~ g 2 T /uj or, if there is a factorized "hard 
thermal loop" involved, even g 2 T 2 /uj 2 . Indeed it is expected that at 0(g 8 ) a contribution 
~ oj5(oj) arises to p e , proportional to the heat capacity and related to the fact that J 
has an overlap with a conserved charge (energy) [25, 26, 24]. At our order, (D(g®), no such 
contributions are present and this is the reason that we have omitted from the outset all 
structures ~ uj n 5(uj) in appendices A and B. 

Even if there are no contributions to G(r) from Lu5(u))-peaks, then how about the infrared 
structures at small but non-zero frequencies, u < gT, discussed in sees. 5.2, 5.3? We note 
that, parametrically, the range uj < m E gives a contribution to eq. (5.27) of magnitude 

fG(r,0) ~ r^LlEM „ g *mlT 2 ~ Otfl*) , (5.32) 

JO TT OJ 

where we made use of the estimate p(uj, 0) ~ g^m^Tu from sec. 5.2. This is formally of higher 
order than the ultraviolet contribution, which is 0(g e T 5 ), but nevertheless more important 
than a full 3-loop analysis, of 0(g s T 5 ). (The estimate in eq. (5.32) also represents the area 
delineated by the transport peak and the w-axis [25].) 
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Figure 6: The ratio (left) and difference (right) of the lattice data from fig. 3(left) of ref. [24] and our 
results from fig. 5(left). The symbols have been slightly displaced for better visibility; the error bars 
have been obtained by varying both results independently within their uncertainties. Lattice errors 
are statistical only; discretization artifacts could affect the difference at the shortest distance, where 
a significant cancellation takes place, the absolute value being about G g /T 5 ~ 50 there. 



After these lengthy qualifications, we apply eq. (5.27) to the results in eqs. (4.1), (4.2), 
(5.11). In accordance with the discussion around eq. (5.32), the significant difference of the 
spectral functions at 0(g 4 ) and 0(g & ) in the infrared regime, visible in fig. 4, converts to a 
rather mild enhancement of the NLO Euclidean correlator around the center of the Euclidean 
time interval in fig. 5. Even though p x is suppressed relative to p e at small frequencies, the 
enhancement is larger in G , because of a contribution from large frequencies. 

Lattice Monte Carlo results for G e /T 5 in pure SU(3) gauge theory have been presented in 
ref. [24]. Multiplying our values by Mac 2 , ^{j^i + j^^) 2 ~ 0.043, we can carry out a direct 
comparison; the results are shown in fig. 6. (We use the 0(g 6 ) rather than the HTL result 
here, because it gives a simpler partial contribution to the transport coefficient, cf. eq. (5.10); 
however the difference is barely visible, cf. fig. 5.) As far as the ratio is concerned, we are quite 
impressed by the good semi-quantitative agreement at short distances. For the difference, the 
good news is that after the subtraction a function is obtained which no longer shows a visible 
short-distance divergence; this implies that a model-independent analytic continuation could 
in principle be attempted. Obviously, however, the significance loss caused by the subtraction 
implies that it remains a tremendous challenge to obtain the per-mille accuracy level that 
might turn out to be necessary for this [6]. (There is the further challenge that, as discussed 
above, a constant r- independent part in G e corresponds to a term ~ u5(u) in the spectral 
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function and does not contribute to the transport coefficient; the resolution has to be good 
enough to show a statistically significant r-dependence in the subtracted correlator.) 

6. Conclusions and outlook 

Even though our study has centered around a specific example, defined in sec. 2 and con- 
cerning the "bulk" and "topology" channels of pure Yang-Mills theory, its wider point, as 
explained in sec. 3, has been to elaborate on a general method for computing next-to-leading 
order thermal spectral functions for composite operators. We have shown how the method 
described yields a relatively simple result, eqs. (4.1)-(4.4), which is well suited for a numerical 
evaluation (fig. 4); as well as how the result can be crosschecked, resummed, and applied in a 
number of ways (sec. 5). Hopefully the step- by-step description of the labour-intensive phase 
of the procedure (appendix A for the most complicated "master" sum-integral appearing at 
2-loop level), will also turn out to lend itself to a number of generalizations; the simplest of 
them may be the inclusion of a non-trivial mass spectrum, which in a specific case has in fact 
already been achieved for many intermediate stages. 

In the course of the computation we have shown that various divergences, sometimes clas- 
sified as being of "soft" , "collinear" , or "thermal" origin, cancel in every of our "master" 
spectral functions separately. That this happens for the sum is in accordance with the gen- 
eral OPE analysis of ref. [20] , showing the absence of infrared divergences up to 4-loop level in 
the regime oj S> vrT, and also with a NLO analysis of the vector current correlator at oj ~ ttT, 
where the cancellation of divergences was verified long ago [15]— [17]. It also agrees with the 
specific discussion of another correlator in ref. [14] for the limit oj S> vrT, which found terms 
proportional to T 2 to be finite at NLO. However, our results are valid down to w ~ irT and 
somewhat below it, and for any single 2-loop master spectral function separately. 

For small frequencies, oj <C ttT, the naive perturbative series does break down. Never- 
theless, as long as g 2 T/ir <C oj <C ttT, it can (probably) still be "repaired" through a Hard 
Thermal Loop resummation (sec. 5.2). The price for the resummation is that in this regime 
only one order of the weak-coupling series is available, and NLO corrections could well be 
substantial [48]. Once oj < g 2 T/ir, even the resummed result loses its validity, and other 
techniques are needed (cf. sec. 5.3); this regime is not addressed in the present study. 

Even if the small-frequency regime of a spectral function suffers from infrared problems, the 
corresponding imaginary-time correlator can still be computed and compared directly with 
lattice simulations, because the contribution emerges dominantly from ultraviolet frequencies, 
w ^ (cf. sec. 5.6). The ultimate use envisaged for such results within thermal QCD is that 
the perturbative imaginary-time correlator can be subtracted from a lattice measurement; 
the remainder could be ultraviolet finite and only sensitive to thermal infrared physics, hence 
showing relatively little "structure"; and that therefore the spectral function corresponding 
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to the remainder could be reconstructed with less model input than without the subtraction. 
That such a program may be feasible is demonstrated by fig. 6(right) for the bulk channel, 
although much further work is needed before the quality of the remainder is satisfactory. 

As promised in the introduction we finally comment on cosmological applications. In the 
present paper the spectral functions were computed at vanishing spatial momentum, i.e. for 
V = (id, 0) in the notation of the introduction, or P = (— ico,0) in the Euclidean notation 
used elsewhere. If, in contrast, the pseudoscalar correlator were computed "on-shell", i.e. for 
P = (—iE p ,p), then it would be directly proportional to the production rate of axions (see, 
e.g., ref. [49]), which are among standard Dark Matter candidates. In the case of a scalar 
density, one could similarly envisage a production rate of dilatons, although in this case it 
might be natural to assume the dilatons to further decay, perhaps into particle-antiparticle 
pairs in a "hidden sector" . Then the dilaton itself could be off-shell, and our result would yield 
the production rate of the hidden-sector pair, in the same way as the spectral function of the 
electromagnetic current in QCD with P = (—uo,0) yields the production rate of a dilepton 
pair with a total energy equal to w (cf. ref. [15]). As far as the corresponding transport 
coefficients are concerned, that related to the pseudoscalar density yields for N c = 2 the 
baryon number violation rate which is certainly important for cosmology [50] (even the case 
= 3 may have some relevance, cf. the discussion in ref. [35] and references therein), 
whereas bulk viscosity could determine e.g. a moduli decay rate [51, 52]. Beyond this, the 
sum-integrals that we have considered are similar to those appearing in other production 
rates, e.g. that of right-handed fermions [53], and so (with the transformation of some lines 
to fermions, which goes simply by changing n B — > — n F ), we envisage that our techniques 
extended to P = (—iE p ,p) might find further applications. In fact, in ref. [14], the results 
were approximated by setting P = (— iM, 0), and then the kinematics is identical to ours. 
We hope to be able to apply our methods to some of these problems in future work. 
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Figure 7: Examples of physical amplitudes contained in eq. (A. 4). Left: "virtual corrections", or 
"factorized integrals". Right: "real corrections", or "phase space integrals". The whole spectral 
function amounts essentially to the absolute value squared of an amplitude, and the virtual process 
contributes through its interference with a tree-level process. (In the case shown, the final answer 
vanishes in the zero-temperature limit, but remains finite for T^O.) 



A. Detailed procedure with the example of p x . 

Let us consider the master sum-integral defined as (this is a slight generalization of eq. (3.15)) 

p6 



HP) - 



Q,R 



Q 2 R 2 [(Q - Rf + A 2 ](Q - Pf{R - Pf ' 



(A.l) 



The corresponding spectral function, evaluated at vanishing spatial momentum, is defined by 



Px ( U )=TmX,(P) 



P->(-i[u>+iO+],0) 



(A.2) 



In the following we give some details on how p x can be given a rapidly convergent 2- 

j 

dimensional integral representation. Since p x is the most complicated structure encountered, 

j 

the same techniques will also apply to the other cases. (The only major simplification in p%. 
is that it is ultraviolet finite so that the analysis can be carried out without ultraviolet reg- 
ularization; nevertheless, the part which contains ultraviolet divergences in other cases will 
already be clearly identified, cf. sec. A.2.) 



A.l. Matsubara sums 

As a first step we carry out the Matsubara sums in eq. (A.l); 1 take the discontinuity leading 
to the spectral function; and set p — > 0. The procedure is similar to that described in some 
detail in refs. [18, 19]. Like in ref. [19] it is useful to "regulate" some of the propagators by 
introducing a small mass parameter which is set to zero in the end (A 2 in eq. (A.l)). To keep 
the situation maximally symmetric, it is convenient to choose the propagator carrying the 
momentum Q — R for this task. Denoting 

E q = q, E r = r, E qr = \/(q - r) 2 + A 2 , (A.3) 

: An alternative strategy might consist of performing first the spatial integrals, which can be done for 
general masses [54] , and by then re-interpreting (some of) the masses as Matsubara frequencies, over which a 
summation is carried out. 
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and not displaying any terms proportional to oj n 5{u) (cf. sec. 5.6), the result from eq. (B.33) 
of ref. [18] can be re-interpreted as 2 



1 



u 6 tt 



q,r 



AqrE, 



qr 



(A.4) 
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The first two structures, with simple 5-function constraints, will technically be referred to as 
"factorized" integrals but correspond physically to virtual corrections; the latter four struc- 
tures, with more complicated (5-constraints, are technically labelled "phase space" integrals 
but correspond physically to real processes (some examples are illustrated in fig. 7). Among 
the factorized integrals a further subdivision into two classes is possible, namely to "pow- 
erlike" integrals in which the (unconstrained) integration proceeds like at zero temperature, 
without any Boltzmann suppression at large momenta (these are virtual vacuum corrections); 
as well as to "exponential" integrals, which contain Bose distributions and are guaranteed to 
be ultraviolet convergent (these are virtual thermal corrections) . We first discuss these three 
classes separately, and then collect together the results. 

It is important to note that the factorized and phase space integrals are both divergent 
(or ill-defined) for A — > 0, because of the appearance of poles in some of the denominators. 



2 We have also rederived the expression from scratch, using the same techniques as in ref. [18] but applied 
directly to the bosonic case, with all propagators taking their tree-level forms. 
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For example in ref. [14] various divergences were classified as being of "soft", "collinear", 
or "thermal infrared" origin (the last coming from at small E). We do not separately 
keep track of the nature of the divergences, but their eventual cancellation, when the real 
and virtual corrections are added together, works out like in ref. [14], and is the reason that 
we can take A — > in the end; however we establish the cancellation for any value of uj/T, 
whereas in ref. [14] only terms of relative magnitude T 2 /u 2 in the OPE-limit were analyzed. 
(Of course, similar cancellations were seen long ago in refs. [15]-[17].) In practice we add as 
a further regulator a principle value prescription for the poles; then changes of integration 
variables, such as q <H> r, are unproblematic and allow to simplify the expressions somewhat. 

A. 2. Virtual vacuum corrections 

Starting with the factorized part, it is immediately clear that the symmetry q o r allows us 
to combine the two first terms of eq. (A. 4). Let us pick the latter case, with 5(uj — 2r) for 
co > 0, as a representative. The integral over r is trivial, 

= (A.5) 

and l/(8r 3 ) from the original "measure" goes over into l/co 3 . In the remaining q- integral we 
substitute the angular variable through E qr : 

Thereby the factorized powerlike integral can be expressed as 

,4 



p 



{q+<% + E qr ){q+<±) ( q -<Z +Eqr )(q-%) 



(A.7) 



The integration range is illustrated in fig. 8. 

Now, we could immediately integrate over the variable E qr in eq. (A.7). It turns out to be 
somewhat simpler to take the limit A — > 0, however, if we change the order of integrations 
(the same change will also be useful in sec. A. 3). This can be achieved through 



roc /V(<7+f) 2 +A 2 ,oo /-f+v^ZA* 

/ dq / dE qr = / dE qr / dq . (A.8) 

Jo V(9-f) 2 + A2 Jx J \%-V E ir-*\ 

Subsequently we rename the inner integration variable to be r, the outer one to be q. Partial 
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Figure 8: The integration range for the factorized integrals, cf. eq. (A. 6), for A = w/10. 



fractioning with respect to the new r, eq. (A. 7) becomes 
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Carrying out the integrals yields 
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where y 7 (^) 2 + A 2 is the value of q for which the lower limit | ^ — \J q 2 — A 2 | of the r-integration 
in eq. (A.9) vanishes. 

Whereas eq. (A. 10) is still exact, we would in the end like to take the limit A — > 0. In 
most of the terms of eq. (A. 10) this can be achieved through a Taylor expansion in X/q, but 
if an integration is sensitive to the range q ~ A then this is not possible, and the integrand 
needs to be kept in its exact form. Analyzing the situation carefully, this can be seen to only 
happen with the very first term, in the vicinity of the lower bound. Thereby we obtain 

(A.ll) 



j 



■5(1 + 2*-) 



s: 



dq 



In 



UJ 

w. 

q + \f<F 



A 2 



q-^qT^tf 



+ 



dq 



In 



(uj + q) uj 



(uj-q)(uj + 2q) 



+ 



'dq 



In 



4{q 2 - uj 2 ) 



Aq 2 



UJ^ 



27 



where indicates that terms of O(AlnA) and smaller have been omitted. Remarkably the 
2nd and 3rd integrals in eq. (A. 11) cancel against each other, so that in total 
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(A.12) 



A. 3. Virtual thermal corrections 

For the factorized exponential integrals we proceed in a similar fashion. The change of 
ordering in eq. (A. 8) and the subsequent renaming of variables is only carried out in terms 
containing the distribution n qr . Thereby the relevant part of eq. (A. 4) becomes 
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The integration over r is again trivial, partial fractioning the latter term like in eq. (A. 9). In 
a few steps we obtain 
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A. 4. Real corrections 



As far as the latter four structures of eq. (A. 4) are concerned, we again rewrite the integration 
measure by substituting the angle between q and r through E qr : 

o pOO pOO pEq r 

— 3 / dq dr dE qr , E% = y/(q± r) 2 + A 2 . (A.15) 

Jo Jo Je- t 



J 



7T 



4qrE qr (4t 



28 



It is also convenient to factor out a common n q n r n qr from the Bose distributions, yielding 

(1 + n gr )(l + n q + n r ) + n q n r = n q n r n qr (e q+r+Eqr - l) , (A. 16) 

n qr (l + n q + n r ) — n q n r = n q n T n qr (e q+r — e Eqr ) , (A. 17) 

n r (l + n q + n qr ) — n q n qr = n q n r n qr (e q+Eqr — e r ) , (A. 18) 

n q (l + n r + n qr ) - n r n qr = n q n r n qr (e r+Eqr - e q ) , (A. 19) 

and to make use of the 5-functions to simplify the denominators. Choosing furthermore 
dimensionless units in the following whereby all variables are expressed in terms of T, and 
noting that, for < A < uj, only four of the ^-constraints get realized, we can rewrite the 
phase space part of the integral as 
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Of the integrals here, the ones labelled (hi) and (iv) are related by the symmetry q <H> r, 
so essentially only three cases remain. We choose to regard the integral over r as an "inner" 
integration, to be carried out first. Then, with some work, the following relations can be 
established for uj > A: 
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Figure 9: Original integration ranges (eqs. (A.21)-(A.23)) for the phase space integrals, for A = w/10. 




co/2 



(i) 



co/2 

(ii) 




Figure 10: Integration ranges after the shifts in eqs. (A.27)-(A.29), for A = w/10. The dashed line 
indicates the diagonal q = r, which can be used to reflect all integrations to the octant q > r. 



where eft is an arbitrary function. The ranges are illustrated in fig. 9. 

It appears advantageous, however, to shift the ranges somewhat, in order to place poles so 
that their eventual cancellation is easier to see. To this end we undertake shifts which turn 
all denominators of eq. (A. 20) into a universal l/{Aqr). The shifts are 

(i) : ?->--g, r -»• - - r , (A.24) 

(ii) : q -»• | + q , r^| + r, (A.25) 

(iii) : q^-^+q, r^^+r, (A.26) 

and the integration ranges of eqs. (A.21)-(A.23) convert into 
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These are illustrated in fig. 10. The expression of eq. (A. 20) becomes 
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2 4q 

Subsequently we can reflect the range q < r to q > r in cases (i) and (ii), as indicated in 
fig. 10; then the integrations only start at q = A/2, and explicit divergences are averted. 

So far we have made no approximations concerning A in the phase space integrals (apart 
from assuming that A < w), but now we again want to send A — > 0. As mentioned this is only 
possible once the factorized and phase space parts are added together, but within each part 
we can extract the "asymptotic" behaviour in this limit. To this end we wish to re-arrange 
the integrations in a form in which a simple divergent part can be computed analytically and 
a more complicated but finite part is left over for numerical evaluation. 

The divergent parts to be subtracted have two separate origins: they are either related 
to a prefactor 1/r, or to the infrared divergent Bose distribution n q + r (or n q _u for the case 
(hi)). To handle n q+r we subtract a term of the form a n q+r / (qr), where a is the "residue" of 
this structure at the origin; the subtracted term can be integrated analytically, by changing 
variables from (g, r) to (x,y) = (q + r, q — r), and carrying out the subsequent integration 
over y. To handle any remaining 1/r-divergences, we subtract a term f3(q)/(qr), where f3(q) 
is the "residue" at r = 0. After these subtractions the remainder integral must remain finite 
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Figure 11: An integration range for which a substitution of the type in eq. (A. 31) can be carried out. 

even in the limit A — > 0. (It is not clear to us whether these subtractions are the unique or 
most elegant ones but they do achieve their goal.) 

As far as the mentioned change of integration variables goes it can be carried out, in a 
range illustrated in fig. 11, as 

f'™d? p»«(«)dr 
1 = / — / —4>{q + r,q-r) 

2 4q 

= / — / h ]<f>(x,v)L _„ , r ,„ \ • (A. 31) 

The upper limit of the y-integration was obtained by eliminating q in favour of x from the 
parametrization of the left boundary of the shaded region in fig. 11: (x,y) = (q+ j^,q~ j^)- 
If the function (j) does not depend on y = q — r, we finally obtain 
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and subsequently rename x to be q. 

Let us illustrate this procedure first with the case (ii) of eq. (A. 30), which is the simplest. 
Symmetrizing in q -H- r, we cancel the 1/2 from the prefactor (the q- integration then starts 
at A/2). From the integrand, +q rim +r n q+r e q+r (l — e w ), we first subtract n q+r together 
with its residue, i.e. 

n 2 u L n q+r (l — e u ) = —(1 + 2n^)n q+r . (A. 33) 

The value of the remainder at r = can be simplified into 

n f +, n % n q e q (l - e") + (1 + 2n« )n q = (1 + 2n f )n q+% . (A.34) 

When both eq. (A. 33) and eq. (A.34) are subtracted, the final remainder simplifies tremen- 
dously, obtaining a form which vanishes for r — > even at small q. In total, this "partitioning" 
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amounts to the identity 



n 



» +q n« +r n q+r e q+r (l-e") = (l + 2n«) 



-n q+r + n q+% -(1 + «,+ «) 



n- 



(A.35) 



The first term of eq. (A. 35) can now be integrated with the help of eq. (A. 32). The second 
term is independent of r so that the principal value integral in eq. (A. 30) evaluates to 
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As mentioned the integral over the third term remains finite when sending A — > 0. Thereby 
the whole expression becomes 
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where "~" is a reminder of having set A — > whenever possible. 

As far as case (iii) goes, the only Bose distribution which diverges at the edges of the 
integration range is n q _^ . Cancelling the factor 2 against 1/2 in the prefactor of eq. (A. 30), 
the residue at q = ^ , r = can be rewritten as 
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When this term is subtracted and r is set to zero, we get 
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Subtracting both terms, we are left over with an integrable remainder; in total, 

(1 + n q - r )(n q - n r+ ± 



n q - r e q -2( e UJ - 1) = (1 + 2n* 



n r n_<£ 

2 



(A.40) 



Now, the principle value integral over the 1/r pole in eq. (A. 30) evaluates to 
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The whole expression then becomes 
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The case (i) is slightly more complicated than (ii) and (hi) because of the non-trivial 
geometry of the integration range (cf. fig. 10); nevertheless, the philosophy remains the same. 
Cancelling the factor 1/2 in eq. (A. 30) for symmetrization, the residue of n q+r at q = r = 
can be written as 
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It will be convenient to also subtract a "vacuum" term in connection with this one, by 
replacing n q+r by 1 + n q+r . When this structure is subtracted and we set r — > 0, we get the 
second subtraction, which can be simplified into — (1 + 2n^)n^ _ q . The final remainder is 
again a function which vanishes as r — > 0. The full partitioning can be expressed as 
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(A.44) 

The first two terms of eq. (A.44), which only depend on q+r, are integrated as in eq. (A. 32). 
The r-integrals over the third term evaluate to 
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Setting again A 



whenever possible, the whole expression becomes 
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A. 5. Collecting the pieces 

Combining now everything together, the first observation is that the factorized powerlike 
term, eq. (A. 12), and the very first term from eq. (A.47), with unity in the numerator, cancel 
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against each other. This implies, in particular, that p x vanishes at T = 0, because these are 
the only terms which are not (necessarily) proportional to a Bose distribution. 

It is more cumbersome to deal with the 1-dimensional integrals containing n q , n\ q _^\, and 
n q+ ". Summing together the terms from eqs. (A. 14), (A. 37), (A. 42) and (A. 47), we obtain 
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(A.48) 



The subsequent step is to take the limit A — > 0. The most subtle in this respect are the cases 
(a) and (b). The various structures appearing there can be simplified as 
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These then lead to 



(a) 
(b) 
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In case (d) we observe that the lower limit of the integration can safely be set to zero in the 
sum. With this knowledge in mind, the left-most terms of (d), (e) and (f) can be summed 



together and the integration variable can be shifted as q — > — ^ + q, to obtain 
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In the right-most terms of (d) and (e) we can shift q — > ^ — q, and in the right-most term 
of (f), we can shift q — > ^ + q, thereby obtaining 
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Summing all of these together with (c), the A's cancel or become harmless, and the regu- 
larization observed in (d), momentarily implicit in eq. (A.54), is rediscovered as a principle 
value integral over l/(q — ^). We thereby obtain the first three rows of eq. (A. 57) below. 

A. 6. Final expression 

Putting together the 2-dimensional integrals from eqs. (A. 37), (A. 42) and (A. 47), and the 
1-dimensional integrals from sec. A. 5, we obtain a convergent expression for p x at A — > 0: 
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For numerical evaluation, it is helpful to combine the first terms of the 2nd and 3rd rows into 
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A. 7. Ultraviolet asymptotics 



As a check of the result obtained, it is useful to verify that it reproduces the known behaviour 
in the limit co S> ttT. From ref. [36], eq. (A. 32), we recall that asymptotically X. reads 
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Setting P = (p n ,0), analytically continuing to Minkowski signature, p n — > —i[co + i0 + ], and 
extracting the spectral function according to eq. (A. 2), yields 
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+ 0(e 2 ), (A.68) 

(A.69) 

+ 0(e 2 ), (A.70) 

(A.71) 

where fx 2 = /2 2 e 7E /47r. Inserting these into eq. (A. 59) we notice that the temperature- 
independent part is of 0(e) and therefore vanishes in the continuum limit, just as we found 
at the beginning of sec. A. 5. The remaining asymptotics reads 

j v ' 32vr 72 \uj 2 ) y ' 

It is easy to check numerically that eq. (A. 57) indeed approaches this form at oj » ttT. 
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B. Main results for the other cases 



Below we list the spectral functions corresponding to the master sum-integrals defined in 
eqs. (3.4)-(3.14). Like in appendix A, all terms containing structures of the type uj n S(uj) are 
omitted from the outset although, as discussed in sec. 5.6, under specific circumstances such 
terms may conspire to give a "zero-mode" contribution; we do not believe this to be the case 
for the observables and perturbative order considered in the present investigation. 



B.l. Pja 

The sum-integral J a is defined as 



J, 



= ^ 



(B.l) 



Since the dependence on p n is polynomial, there is no cut, and the corresponding spectral 
function vanishes: 

Pj(u) = 0. (B.2) 



B.2. 



The sum-integral J h is defined as 



J h = 
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Q 2 {Q-P? 



(B.3) 



Carrying out the Matsubara sum, taking the cut, and setting subsequently p — > 0, we get 
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Making use of the usual dimensionally regularized integration measure, the integral over the 
Dirac-5 yields 



tt5(uj — 2q) = 
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l + e(ln^ + 2) +0(e 



and the complete expression reads (for uj > 0) 

w 4 -2e 

p 7 (oj) = — *1 (i + 2n-) 

' J * y ' 16vr v 2 ' 



I : In ^ • 2 +0{c 2 ) 



(B.5) 
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Thermal effects are seen to be exponentially suppressed here, so we nicely match the asymp- 
totic OPE-result as determined in ref. [36], 
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B.3. p Xa 

The sum-integral X a is denned as 



Since there is no dependence on P, the corresponding spectral function vanishes: 

0z» = O. (B.9) 

B.4. PXb 

The sum-integral X b is defined as 

lh = t QR Q 2 R 2 {R-P? ■ (R10) 
This factorizes into the same form as J h in eq. (B.3), multiplied by the well-known 



L-L 2 



— = — + 0(e). (B.11) 



Within the 2- loop contributions we only need to go to 0(e°), so the result reads 

to > ) = "lS 1 + 2n t)- (B.12) 
For uj » 7rT this agrees with the asymptotic OPE- form from ref. [36], 



, .2rp2 /rpQ 



B.5. p Xc 

The sum-integral X c is defined as 



- 1^' 



Since the dependence on p n is polynomial, there is no cut, and the corresponding spectral 
function vanishes: 

fe» = 0. (B.15) 
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B.6. p T 

The sum-integral I d is denned as 

Xd = t QR Q 2 R\R-P? • (R16) 

This factorizes into the tadpole in eq. (B.ll), times a 1-loop sum-integral. For the latter, a 
possible trick is to write 
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The result for the Matsubara sum here (after taking the cut and setting p — > 0) can be read 
directly from eq. (B.4), 
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where E 2 = r 2 + A 2 and the spatial integral was performed for oj > 0. Taking the derivative 
in eq. (B.17) and setting A — > finally yields 

uj 2 T 2 2uj 2 f°° 

fe d H = T92^( 1 + 2n t) = (4^( 1 + 2n f) J, dqn « q > (R19) 

where we have also given an integral representation, in order to facilitate comparison with 
the other spectral functions below. The result matches the OPE-behaviour from ref. [36], 

2 T 2 /T 6 \ 
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192vr \uj\ 
B.7. Plc 

The sum-integral Z c is defined as 

Xc 55 t QR Q 2 R 2 {Q-R?- (R21) 

Again the result is a polynomial in p n and has no cut; in fact even before the cut this massless 
sum-integral vanishes exactly in dimensional regularization in any dimension [45]. Therefore, 
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B.8. p Xf 

The sum-integral X f is denned as 



2f = 



Q,R 



Q 2 [(Q-Ry + \*](R-P) 2 



(B.23) 



After the Matsubara sums the expression for the corresponding spectral function reads 
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Evidently this sum-integral only contains a phase-space part, and the expression correspond- 
ing to eq. (A. 30) becomes 
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Subsequently we can reflect the range q < r to q > r in cases (i) and (ii). The subtractions of 
the potentially divergent parts are carried out like in eqs. (A. 35), (A. 40), (A.44); the integrals 
over y and r are trivial in these terms. 

Once the pieces are collected together, we can again identify a "vacuum part", a "1- 
dimensional integral" , and a "2-dimensional integral" . The contribution to the vacuum part 
from the phase space integrals reads 
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This yields the first row of eq. (B.27). The 1-dimensional integrals have no delicate divergences 
here. Setting A to zero, we obtain 
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10 Jo n r 

For uj » 7rT it can be checked numerically that this goes over to the known [36] ultraviolet 
asymptotics 

(B.28) 
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B.9. p T 

The sum-integral l g is defined as 
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(B.29) 



After the Matsubara sums the expression for the corresponding spectral function reads 
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The two terms can be combined through q <H> r if we first introduce a principle value pre- 
scription. Choosing the first one as a representative, the r-integral can be carried out like in 
eq. (B.5), and we obtain 
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The thermal part (with n q ) is finite as it stands, whereas the vacuum part may be transformed 
into a well-known four-dimensional integral: 
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This can be verified by carrying out f ^ on the right-hand side. In total, omitting 0(e), 
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For oj » 7rT this approaches the OPE-result [36 
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B.10. /9 X 



The sum-integral X h is defined as 
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After the Matsubara sums the expression for the corresponding spectral function reads 
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The main difference to the procedure described in appendix A is the handling of the 
"factorized powerlike" integrals, which are ultraviolet divergent, like in eq. (B.32). The 
expression reads 
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The r-integral can be evaluated like in eq. (B.5), whereas the q-integral can be re-expressed 
as a four-dimensional vacuum integral: 
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This can be verified by carrying out J ^ on the right-hand side. The value of the integral is 
familiar; in fact, since 0(D) rotational invariance implies that the result depends on A 2 and 
R 2 only, and R 2 = according to eq. (B.38), the simplest way to derive it is just by setting 
R = 0. In any case, 
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The "factorized exponential" integrals can be worked out like before; the expression corre- 
sponding to eq. (A. 14) becomes 
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For the phase space integrals the expression corresponding to eq. (A. 30) reads 
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Subsequently we can reflect the range q < r to q > r in cases (i) and (ii) . The subtractions of 
the potentially divergent parts are carried out like in eqs. (A. 35), (A.40), (A.44); the integrals 
over y and r are a bit more complicated than with Zj but can be carried out. 

Once the pieces are collected together, we can again identify a "vacuum part", a "1- 
dimensional integral" , and a "2-dimensional integral" . The contribution to the vacuum part 
from the phase space integrals reads 
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Summing together with eq. (B.40), we obtain the first row of eq. (B.44). The most tedious 
work concerns the handling of the 1-dimensional integrals; however in principle everything 
proceeds like in sec. A. 5, and eventually the A's cancel or can safely be put to zero. Altogether 
we obtain 
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For oj 3> 7rT it can be checked numerically that this goes over to the known ultraviolet 
asymptotics [36] 

B.ll. p x 

The sum-integral X, is defined as 
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At zero spatial momentum, 4(Q ■ P) 2 = ^q\p 2 n = 4(Q 2 — q 2 )p 2 n , so that 
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The Matsubara sums here are equivalent to those in I h and I h , and we obtain 
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(q + r - E qr )(q-r - E qr ) 
n r (l + n q + n qr ) 

(q - r + E qr )(q + r + E qr ) 
n q (l + n r + n qr ) - n r n qi 



E qr )(q + r-E, 
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Like with p x the "factorized power like" integrals are ultraviolet divergent; the expression 

h 



reads 



u 2 qir 



4r 2 E, 



S(u-2r)[ — + 



qr 



+ 



(1 + 2n r ) . (B.49) 



q q + r + E qr -r + q + E qr t 

The r-integral is the same as in eq. (B.5), whereas the non-trivial q-integral can be re- 
expressed as a four-dimensional vacuum integral: 



+ 



+ q + E qr -% + q + E, 



qr 



Q Qi[(Q- R y + \2] 



(B.50) 



i?=(|i,|e r ) 



This can be verified by carrying out J ^ on the right-hand side of eq. (B.50). Once expressed 
this way, we can make use of 0(D) rotational covariance, in order to reduce the "tensor" 
structure in the numerator to scalar integrals. The result is then a function of A 2 and R 2 
only; given that R 2 = 0, we need the leading term in a Taylor-expansion around this point. 
A few steps lead to 



Ic 



Q Q2[(Q-Rf + X 2 } 



U 2 ^ 
12(4vr) 2 



i+^ + i+oie: 
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and in total 



12(4vr) 3 



7 



(l + 2n f )(± + ln^ + ln^ + ^+0( e ) 



(B.52) 



The "factorized exponential" integrals can be worked out like before; the expression corre- 
sponding to eq. (A. 14) becomes 



(47T) 



2gw + A 2 



2gw - A 2 



In 



•oo 

- / dqn q x 

,2 



+ ^-f) 2ln 



A 2 



For the phase space integrals the expression corresponding to eq. (A. 30) reads 



(4vr) : 



■(e--l) 



(i) - 



dg 



2u 



dr F 



(?-f) 2 



uj\2-i 



f 00 r°° Ua + -) 

(h) +y d,y^drp ii± 2 - 



n| +g n- +r n g+r e 



(hi) + / dq 



j(2q-u)-\' ! 



4q 



dr F 



(B.53) 
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Subsequently we can reflect the range q < r to q > r in cases (i) and (ii). The subtractions of 
the potentially divergent parts are carried out like in eqs. (A. 35), (A. 40), (A. 44); the integrals 
over y and r are a bit more complicated than with I j but can be carried out. 

Once the pieces are collected together, we can again identify a "vacuum part", a "1- 
dimensional integral" , and a "2-dimensional integral" . The contribution to the vacuum part 
from the phase space integrals reads 



^(l + 2n f )/ A 2 d,[(,-|) 2 ln 



q-^fqT^)2 



uq + - q' 



^ H o \(, A 2 13 \ 

— ; — r^-fl + 2n<A I In — ^ H 

12(4vr) 3 1 2 ' \ co 2 6 J 
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Summing together with eq. (B.52), we obtain the hrst row of eq. (B.56). The most tedious 
work concerns the handling of the 1-dimensional integrals; however in principle everything 
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proceeds like in sec. A. 5, and eventually the A's cancel or can safely be put to zero. Altogether 
we obtain 



(47t) 3 Px.M 
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2cj(1 + 2n<±) 
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Numerically it can checked that for u » ttT this approaches the known asymptotics [36] 
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C. Hard Thermal Loop resummation 

We give here some details for the Hard Thermal Loop (HTL) resummation discussed in 
sec. 5.2, which is important when oj ~ gT and could allow us to postpone the breakdown of 
the perturbative series down to lower frequencies (oj < g 2 T/i:). 

In general, HTL resummation is needed both for propagators and vertices; for instance, the 
coupling of a photon to quarks and gluons gets modified by a correction oc g 2 T 2 and needs 
to be taken into account at leading order if the photons are soft [43]. On the other hand, 
the coupling of a heavy quark to gluons does not get modified even at NLO within HTL 
perturbation theory [48]. For us relevant is the coupling of a dilaton or an axion to gluons. 
For "hard" axions it can trivially be argued that no vertex modification is needed [55] but here 
we are concerned with "soft" dilatons or axions (oj ~ gT, p = 0). It is nevertheless our belief 
that no vertex modification is needed, neither for dilatons nor for axions. Concrete indications 
in this direction are that even without a vertex correction we find results which exactly cancel 
the leading infrared divergences from the QCD results in both channels (cf. sec. 5.2), and 
that the results satisfy the expected sum rules (cf. last paragraph of sec. 5.4). In addition, 
we have checked the absence of vertex corrections through explicit 1-loop computations. 3 
We assume that simple theoretical arguments could also be given for the absence of vertex 
corrections but at the time of writing these evade us. No further details are provided here 
because, strictly speaking, it would not even be necessary for us to treat the infrared regime 
correctly. 



C.l. Contractions 



With the logic outlined above, the only ingredient needed from the HTL theory is the gauge 
field propagator, which takes the form 



(A°{X) A b u (Y)) = 5 ab ^ e iQ < x ~ Y) 



+ 



+ 



Q 2 + U T (Q) Q 2 + U E (Q) Q 4 
where £ is a gauge parameter. The two independent projection operators are defined by 

Q fiQv 



(C.l) 



f^(Q) = [sa - ^) , pJUq) = V 



Q 2 



KAQ) , 



whereas the self-energies read 
U T (Q) 



mt 



D 



Ik 



Q 2 



qz 



q iq n + qz 



n £ (0) = 



qz 



iq n + qz 



(C.2) 



(C.3) 



3 According to ref. [56], considering a three-leg process in Feynman gauge Hard Thermal Loops could 
arise from structures of the types % K k-i^-p^k-q-)* > K^ic-p)^ ' as weu as $kt?i- I n tne dilaton case 
individual graphs do give contributions of the latter two types but these cancel in the sum; in the axion case 
we find no such contributions. 
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with the Debye mass parameter given by 



m 



2 =g 2 N c (D-2) 2 Jj^l 



(C.4) 



For D = 4 this reproduces the expression given in eq. (5.13). 

In the main body of the text, a next-to-leading order computation was described. For 
uj ~ gT, however, the series representation breaks down. Through the HTL method we then 
resum infinite orders to a single contribution, which is of the leading order within the HTL 
power counting. At this order, all the integrals met are ultraviolet finite, so that we could 
set D = 4; for generality, however, we mostly show D in its original form. 

When the propagator takes the form in eq. (C.l), the contractions leading to eqs. (3.2), 
(3.3) need to be carried out anew. Because of the non-trivial structure of eq. (C.l), this is 
somewhat tedious even at leading order. The problem gets simplified if we set p — > from 
the outset; then 



Also, a combination appearing frequently, particularly in the x-channel, reduces to Q 2 P 2 — 
(Q • Pf = q 2 pl. Thereby we obtain (making also use of [Q ■ (Q - P)] 2 = Q 2 (Q - P) 2 - 
[Q 2 P 2 — (Q ■ P) 2 } in order to identify a common part in both channels) 



where we have defined A T (Q) = Q 2 + Il T (Q) and A E (Q) = Q 2 + U E (Q). 
C.2. Naive HTL computation 

As discussed around eq. (5.11), we need both a "naive" and a "resummed" version of the 
HTL result. We start with the naive computation; this means that va\ is treated as a small 
quantity, of 0(g 2 ), and the HTL-result is "re-expanded" in g 2 up to NLO. 

4 We have employed a notation valid in an arbitrary dimension, with the integration measures defined as 
(d = D - 1 = 3 - 2e) 



P^lAQ) = Q^UQ) = • 



(C.6) 





(C.8) 




(C.5) 



such that f = f f and f 1 = 1. 
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Now, since the HTL theory (correctly) modifies physics only at soft scales, one should be 
careful not to apply HTL resummation to parts of the phase space where it is not needed. 
In real-time computations, it is notoriously difficult to identify the relevant parts in advance. 
Therefore, in the following we carry out the HTL resummation on all scales, whereby the 
HTL results show dependence on two scales, both the "soft" m E as well as the "hard" nT. 
The unwanted dependence on irT is easily dropped a posteriori, cf. text below eq. (C.25). 

With this philosophy, we start with the pseudoscalar (%) channel, which is simpler (cf. 
eq. (C.8)). It will be checked afterwards that the additional terms in eq. (C.7) give no 
contribution in the "naive" limit (apart from "contact terms" of the type discussed in sec. 5.4). 

First of all, at leading order, A T (Q) = Q 2 . Then the structure in eq. (C.8) is almost the 
same as the master sum-integral J h of eq. (B.3), and we get 



Im 



= [ 

Q 2 {Q - P) 2 J P-».(-i[ w +iO+],0) Jq 4 



5(u - 2q) - 5{uj + 2q) (1 + 2n q ) 



Adding the factor 2{D — 2)g A and setting D->4we get 

„HTL /, ,\ 4 4 

^|)^ (1 + 2 „ f)+0(9 «), ( C,0) 

which indeed agrees with the leading-order term in eq. (4.2). 

At NLO, we can expand A T to first order in n T . Substituting subsequently Q — ^ P — Q 
yields 



G» TL (P) 



-16<U4(£>-3) 



«"-*>A4 (en) 



d(<? 6 ) " Q 

where Ii T is to be inserted from eq. (C.3). There are two structures; for the simpler one we 
proceed in analogy with eqs. (B.17), (B.18): 



■2 



T Q Q 4 (Q - P) 2 2 &\ 2 f Q [Q 2 + X 2 ][(Q - P) 2 + X 2 } 



2 dA 2 



J ^S( U -2E q ){l + 2n Eq ) 



1 _d_ f K 



q Q 

2 A \2 



(l + 2nu)0(u-2\) + O(e) 



2 d\ 2 { 64-koj 

(l + 2n«) . (C.12) 



4(4vr) 
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The second term is a bit more complicated, and it is useful to introduce an infrared regulator: 



't Q I z Q 2 (Q- 

- -L 

- -I I 

Jz,q JO 



P) 2 iq n + qz 



(q n - Pn) 2 + q 2 ql + E 2 iq n + qz 



(36., 



ffh ~\~Pn <?n 



1 



3 

dr 
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qz 



( r n + q 2 ){ql + E T) 2 \iq n + qz -iq n + qz 



2 ^2 ^i'pn 



tr n r 



E 2 - q 2 z 2 



^ r 2 +q 2 ^ 



-IQnT 



ll + q 2 z 2 q 2 n + E 2 



(C.13) 



where E 2 = q + A ; we symmetrized in z ■<->■ — z; represented the Kronecker-<5 as an integral; 
and partial- fractioned the dependence on q 2 . The sums can be carried out, e.g. 



n n 



r„ 



r 2 + q 2 2q 



e qT + e q 



3-r) 



(C.14) 



and subsequently also the integral over r, yielding 



I(Pn) = 



q 2 z 2 



n„ 



z,q 

n qz 
qz 



X 2 + q 2 (l-z 2 )4q 



+ 



e q(l+z)l3 _ ^ 



1 



+ 



1 



-iPn + q(l + z) ip n + q(l + z\ 

1 1 

+ 



Eq [\ J \~iPn + 

+ (> - e E ' 



ip n + q(l - z) ip n + q(l - z) 

1 



+ 



q + E g ip n + q + E q 

1 1 

+ 



(C.15) 



. -iPn + q~ E q ip n + q - Eg / 

We can now set p n —> + i0 + ] and take the imaginary part; for oj > A > 0, only three 
channels contribute. If we also substitute z —> —z in one of them, and make use of the 
5-constraints in order to re-express the arguments of the Bose distributions, the expression 
reduces to 

■(l + Tig + Hoj-g) 



lml(-i[u} + i0 + ]) = -j f 

J Z. 



z 



:,q \ 2 +q 2 {l-Z 2 ) 

- 5(uj — q — qz) — J- 8(uj — q — E q ) 



L z 



E n 



The integrals over z can be carried out: 



1 f +1 , 2z5{uj-q-qz) 
X 2 + q 2 (l-z 2 ) 



dz 



2 7-1 X 2 + q 2 {l-z 2 ) 



^-q)e(q-%) 
q 2 (X 2 + 2u;q-oj 2 ) ' 

2g £7, + q 



(C.16) 

(C.17) 
(C.18) 
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The q- integral over 5(u) — q — E q ) yields 



f°° E 
/ dq<f>(q,E q )5(uj-q-E q ) = — ^<t>(q,E q 

JO Q + E q 



1- 9„, ■ ^ 



(C.19) 



whereas the other (/-integral is split into the ranges , cu) and (uj, oo). In the former range, we 
add and subtract 1 + 2n^ from the Bose distributions, whereas in the latter range, we write 
1 + — q — Tlq — ijj in order to have a positive argument. In total, setting A — > whenever 
possible, 



Im/(-i[w + 



1 

'8tt 



, uj — q , n. 
dq ^Jq—^^ + n ^- 2n V 



(C.20) 
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uj — q , \ 

' 1 ^2^17^-^ 



+ / dq 



E, 



A 2 + 2w<? - w 2 v 2 ' 



A. + 1 i n -^g ~g 



2cj 



The last two terms contain a logarithmic divergence, but their sum is finite: 

■ 2 1 



s: 



dq 



uj — q 



A 2 + 2ujq - uj 2 



E n 



q + E q \E q 2 E q + q 
If we also write 

n q + n^-g - 2n^ 



(l + 2n f ) 



n g (l + n. 



n n -n 



q—U) 



(l + 2n f ) 



+ n. 



q—u) J 



nz 



(C.21) 
(C.22) 

(C.23) 
(C.24) 



and combine then with the contribution in eq. (C.12) as well as the leading-order term from 
cq. (C.10), the full result becomes 



(47T) 



■(l + 2n^)\uj 4 - iuj 



I _ 2 ("tig UJ ~ q re 9( 1+n ^-g) 
4 7^ 2wg — oj 2 



Ju 2wg - w 2 



Vf 

+ 0(g 8 ). (C.25) 



In a systematic HTL-computation, we need furthermore to replace the distribution functions 
through their "classical" limits, (l+)n q — > T/q etc.; then the integrals are trivially doable 
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and we get eq. (5.12). (For notational simplicity we however always display the prefactor 
1 + 2m± in its "quantum" form.) 

We end by showing that additional terms in Gg TL , on the second row of eq. (C.7), give no 
cu-dependent contribution in the "naive" limit. At leading order this is clear: both A E (Q) and 
A T (Q) get replaced with Q 2 , so that we are faced with 2g 4 $ Q l, which vanishes in dimensional 
regularization. At NLO, expanding in II^, II T and substituting Q — > P — Q, we obtain 

grCg) , G^(P) r r u T (Q) 

Ad A c 2 e ^ 16d A c 2 (D - 3) T \ Q 2 ' Q 2 



t h • (a26> 



where we inserted the self-energies from eq. (C.3). This amounts to a "contact contribution": 
as discussed around eqs. (5.20), (5.21), at 0(<? 6 ) the two channels do differ in this respect. 
However, since eq. (C.26) is P-independent, there is no contribution to the spectral function. 

C.3. Resummed HTL computation 

When eqs. (C.7), (C.8) are kept in their unexpanded forms, the resulting spectral functions 
can be expressed as convolutions of two "elementary" spectral functions. More precisely, if 
we make use of the spectral representation, 

I _ f°° dwi p T (ui,q) I _ [°° du 2 P T (u2,q) (C 27) 



& T (Q) y_oo vr wi-iq n ' A T (Q - P) 7_ 00 vr u 2 - i{q n - p n ) 
then (cf. e.g. ref. [43]) 

1 



1 ^ A T (Q)A T (Q - P) j PM _ i[u}+iQ+l0) 

= (l + 2n f ) r^p T (q°, q) p T (u - q°,q) . ( C .28) 

2 

Let us recall a derivation. Inserting eq. (C.27) into the left-hand side of eq. (C.28), the sum 
can be carried out like in eq. (C.13): 

4^ [u>i - iq n ][u2 ~ i(q n ~ Pn)} ^ ^ (wi - iq n )(^2 - ir n ) 



r/3 p -iqnT 

/ dre'^rV- TV 

Jo „ ^i - % q n „ 



HTt 1 n 
V v v ' 

= n wl n W2 ■ . (C.29) 
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The summations over q n and r n here are "marginally" convergent; their results are dis- 
continuous at r = and r = /3, but there are no Dirac-5's, so the subsequent integra- 
tion is unproblematic. Setting p n — > —i[uj + i0 + ] and taking the imaginary part yields 
Im{l/[o; + LU2 — u>\ + i0 + ]} = — tt5(uj + uji — so that the integral over coi can be carried 
out. Afterwards we write 

n u)+uj2 n uJ2 ^ 2 ^ — = n^-i-^ (l ~l~ n ui2^) -0 

+ 1 

= -n^+^n-^ . (C.30) 

2 

Substituting W2 — > —q° and making use of Px(—q ,q) = —Pt(q i<i) we recover eq. (C.28). 

In order to make use of eq. (C.28), the objects appearing should be regular enough to 
allow for a faithful spectral representation; as discussed in sec. 5.4, this requires that the 
corresponding Euclidean correlators vanish as p n — > oo. This is not the case with the terms 
on the second row of eq. (C.7); however, adding and subtracting ITs in the numerator, and 
shifting Q — >■ P — Q where appropriate, it is easy to see that apart from P-independent terms 
the result can be rewritten as 



2(D - 2)gV n f 



■2 



4d A <% v JqA t (Q)A t (Q - P) 

+ 29 T Q \A E (Q)A E (Q-P) +iD 2) A T (Q)A T (Q-P)/ +C ° nSt - (C ' 31) 

Now all individual "propagators" decrease as a function of the four-momentum. 

From eqs. (C.8), (C.27), (C.28), (C.31) we observe that the following spectral functions are 
needed: 



l { A T (q n ,q) 



p T ( q °,q) = Im\ A }_ ^ \ (C.32) 



Pt^Q)^\^M:\ , (C33) 

p E tf,q)=1m\%^\ . (C.34) 

To determine these we recall that, for D — > 4, the z-integrals in eq. (C.3) can be carried out: 



U E (-i[q +iO + ],d) = m\ 



o\2 



q 2 



q° </» • q I /I)' , 
1 "2^ ln g 0- g + i0+ ^ (C - 36) 
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Denoting 

a straightforward computation leads to 

r T (r?) 



MQ) = { 



Pt(Q) 



r}= — 
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H < i , 



p B (Q) 



S2(Q)+r2( r/ ) ' 
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Ifjl > 1 , 



where we have introduced the well-known functions [57]-[59] 
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For integrals over the poles it is useful to note that 
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The approximate pole positions can also be worked out in specific limits: S T (Q) 
V ~ T7& + IS for « < m E and 7/ « 1 + $ for 9 » m B , whereas S £ (Q) = 



^ 5m E 
3^9 



77 « -EE- + 2g for q <C m E and r/ « 1 + 2e" 2 ^ / m E+!) for g > m B . 
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Figure 12: An illustration of the phase space relevant for HTL-resummation. The shaded areas and 
the dash-dotted blue lines indicate regions in which at least one of the spectral functions has non-zero 
support; the final contribution emerges from the part where both are simultaneously non-zero, which 
happens within the hashed areas. For oj > 2ojp\, where uj p i = m E /\/3, this leads to pole-pole ("p-p"), 
pole-cut ("p-c") as well as cut-cut ("c-c") contributions (cf. ref. [43]). As indicated by the arrow, 
symmetry allows us to restrict the integration to the range q° > lu/2. 



With the spectral functions at hand, we are left to insert them to the imaginary parts of 
eqs. (C.8), (C.31), rewritten through eq. (C.28): 



7T5f 4 (l + 2n|) r 16 



4d A d 
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(C.49) 
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.0 „ („Q „\„ (,, „0 ^^ 9° u-q° 



p t (q ,q)pt( U} ~ Q 



11 z 



(C.50) 



where we have identified the same prefactor as in eq. (C.25). Given the structures of 
eqs. (C.38)-(C.40), these integrals include three types of contributions, illustrated in fig. 12 
(the situation is similar to that in ref. [43]). Taking p^ TL as an example, the "pole-pole" term, 
which only contributes for oj > 2a> p i, where uj p \ = m E /\/3 denotes the "plasmon frequency", 
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amounts to 
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(C.51) 



,q)=0 



Here we made use of the fact that along the curve on which S T (Q) = 0, ^ = — f g ^ T 

For the "pole-cut" contribution, making use of the reflection symmetry with respect to the 
— - (cf. fig. 12), we similarly get 
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(C.52) 



S T (gO,«j)=0 



where the value of g m i n is determined from T, t (lj — g m in, (/min) = for u > w p i, and S T (o; + 
Qmm,Qmm) = for w < w p i. Finally, the "cut-cut" term reads 



7T 2 ./<£ ./<£ 



r T (g°,g) 



r T (^-g°,g) 



S 2 (g°, g) + r^gO, g) ( w - g°, g) + r^(oj - g°, g) n? 



(C.53) 



In a systematic HTL-computation, we are furthermore to replace the Bose distributions in 
cqs. (C.52), (C.53) through their "classical" limits, n q0 — > T/q° etc., so that n q0 n uj _ q0 /rvl — > 

V(oj-q u ) ' Thereby the HTL-results only feel the scale ^ ; this is in analogy with the thermal 
corrections, which only depend on 
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